home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 11 / Cream of the Crop 11-1.iso / math / ast51src.zip / CHARTS3.C < prev    next >
C/C++ Source or Header  |  1995-12-31  |  27KB  |  720 lines
  1. /*
  2. ** Astrolog (Version 5.10) File: charts3.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1995 by Walter D. Pullen
  6. ** (Astara@msn.com). Permission is granted to freely use and
  7. ** distribute these routines provided one doesn't sell, restrict, or
  8. ** profit from them in any way. Modification is allowed provided these
  9. ** notices remain with any altered or edited versions of the program.
  10. **
  11. ** The main planetary calculation routines used in this program have
  12. ** been Copyrighted and the core of this program is basically a
  13. ** conversion to C of the routines created by James Neely as listed in
  14. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  15. ** available from Matrix Software. The copyright gives us permission to
  16. ** use the routines for personal use but not to sell them or profit from
  17. ** them in any way.
  18. **
  19. ** The PostScript code within the core graphics routines are programmed
  20. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  21. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  22. **
  23. ** The extended accurate ephemeris databases and formulas are from the
  24. ** calculation routines in the program "Placalc" and are programmed and
  25. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  26. ** (alois@azur.ch). The use of that source code is subject to
  27. ** regulations made by Astrodienst Zurich, and the code is not in the
  28. ** public domain. This copyright notice must not be changed or removed
  29. ** by any user of this program.
  30. **
  31. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  32. ** X Window graphics initially programmed 10/23-29/1991.
  33. ** PostScript graphics initially programmed 11/29-30/1992.
  34. ** Last code change made 12/27/1995.
  35. */
  36.  
  37. #include "astrolog.h"
  38.  
  39.  
  40. /*
  41. ******************************************************************************
  42. ** Multiple Chart Scanning Routines.
  43. ******************************************************************************
  44. */
  45.  
  46. /* Search through a day, and print out the times of exact aspects among the  */
  47. /* planets during that day, as specified with the -d switch, as well as the  */
  48. /* times when a planet changes sign or direction. To do this, we cast charts */
  49. /* for the beginning and end of the day, or a part of a day, and do a linear */
  50. /* equation check to see if anything exciting happens during the interval.   */
  51. /* (This is probably the single most complicated procedure in the program.)  */
  52.  
  53. void ChartInDaySearch(fProg)
  54. bool fProg;
  55. {
  56.   char sz[cchSzDef];
  57.   int source[MAXINDAY], aspect[MAXINDAY], dest[MAXINDAY],
  58.     sign1[MAXINDAY], sign2[MAXINDAY], D1, D2, occurcount, division, div,
  59.     fYear, yea0, yea1, yea2, i, j, k, l, s1, s2;
  60.   real time[MAXINDAY], divsiz, d1, d2, e1, e2, f1, f2, g;
  61.   CI ciT;
  62.  
  63.   /* If parameter 'fProg' is set, look for changes in a progressed chart. */
  64.  
  65.   ciT = ciTwin;
  66.   fYear = us.fInDayMonth && (Mon2 == 0);
  67.   division = (fYear || fProg) ? 1 : us.nDivision;
  68.   divsiz = 24.0 / (real)division*60.0;
  69.  
  70.   /* If -dY in effect, then search through a range of years. */
  71.  
  72.   yea1 = fProg ? Yea2 : Yea;
  73.   yea2 = fYear ? (yea1 + us.nEphemYears - 1) : yea1;
  74.   for (yea0 = yea1; yea0 <= yea2; yea0++) {
  75.  
  76.   /* If -dm in effect, then search through the whole month, day by day. */
  77.  
  78.   if (us.fInDayMonth) {
  79.     D1 = 1;
  80.     if (fYear) {
  81.       Mon2 = 1; D2 = DayInYearHi(yea0);
  82.     } else
  83.       D2 = DayInMonth(fProg ? Mon2 : Mon, yea0);
  84.   } else
  85.     D1 = D2 = Day;
  86.  
  87.   /* Start searching the day or days in question for exciting stuff. */
  88.  
  89.   for (Day2 = D1; Day2 <= D2; Day2 = AddDay(Mon, Day2, yea0, 1)) {
  90.     occurcount = 0;
  91.  
  92.     /* Cast chart for beginning of day and store it for future use. */
  93.  
  94.     SetCI(ciCore, fYear ? Mon2 : Mon, Day2, yea0, 0.0, Dst, Zon, Lon, Lat);
  95.     if (us.fProgress = fProg) {
  96.       is.JDp = MdytszToJulian(Mon2, DD, yea0, 0.0, Dst, Zon);
  97.       ciCore = ciMain;
  98.     }
  99.     CastChart(fTrue);
  100.     for (i = 1; i <= cSign; i++) {
  101.       cp2.cusp[i] = house[i];
  102.       cp2.house[i] = inhouse[i];
  103.     }
  104.     for (i = 1; i <= cObj; i++) {
  105.       cp2.obj[i] = planet[i];
  106.       cp2.dir[i] = ret[i];
  107.     }
  108.  
  109.     /* Now divide the day into segments and search each segment in turn. */
  110.     /* More segments is slower, but has slightly better time accuracy.   */
  111.  
  112.     for (div = 1; div <= division; div++) {
  113.  
  114.       /* Cast the chart for the ending time of the present segment. The   */
  115.       /* beginning time chart is copied from the previous end time chart. */
  116.  
  117.       SetCI(ciCore, fYear ? Mon2 : Mon, Day2, yea0,
  118.         DegToDec(24.0*(real)div/(real)division), Dst, Zon, Lon, Lat);
  119.       if (fProg) {
  120.         is.JDp = MdytszToJulian(Mon2, DD+1, yea0, 0.0, Dst, Zon);
  121.         ciCore = ciMain;
  122.       }
  123.       CastChart(fTrue);
  124.       for (i = 1; i <= cSign; i++) {
  125.         cp1.cusp[i] = cp2.cusp[i]; cp1.house[i] = cp2.house[i];
  126.         cp2.cusp[i] = house[i];  cp2.house[i] = inhouse[i];
  127.       }
  128.       for (i = 1; i <= cObj; i++) {
  129.         cp1.obj[i] = cp2.obj[i]; cp1.dir[i] = cp2.dir[i];
  130.         cp2.obj[i] = planet[i];  cp2.dir[i] = ret[i];
  131.       }
  132.  
  133.       /* Now search through the present segment for anything exciting. */
  134.  
  135.       for (i = 1; i <= cObj; i++) if (!ignore[i] && (fProg || FThing(i))) {
  136.         s1 = SFromZ(cp1.obj[i])-1;
  137.         s2 = SFromZ(cp2.obj[i])-1;
  138.  
  139.         /* Does the current planet change into the next or previous sign? */
  140.  
  141.         if (!ignore[i] && s1 != s2 && !ignore[0]) {
  142.           source[occurcount] = i;
  143.           aspect[occurcount] = aSig;
  144.           dest[occurcount] = s2+1;
  145.           time[occurcount] = MinDistance(cp1.obj[i],
  146.             (real)(cp1.dir[i] >= 0.0 ? s2 : s1) * 30.0) /
  147.             MinDistance(cp1.obj[i], cp2.obj[i])*divsiz + (real)(div-1)*divsiz;
  148.           sign1[occurcount] = sign2[occurcount] = s1+1;
  149.           occurcount++;
  150.         }
  151.  
  152.         /* Does the current planet go retrograde or direct? */
  153.  
  154.         if (!ignore[i] && (cp1.dir[i] < 0.0) != (cp2.dir[i] < 0.0) &&
  155.           !ignore2[0]) {
  156.           source[occurcount] = i;
  157.           aspect[occurcount] = aDir;
  158.           dest[occurcount] = cp2.dir[i] < 0.0;
  159.           time[occurcount] = RAbs(cp1.dir[i])/(RAbs(cp1.dir[i])+
  160.             RAbs(cp2.dir[i]))*divsiz + (real)(div-1)*divsiz;
  161.           sign1[occurcount] = sign2[occurcount] = s1+1;
  162.           occurcount++;
  163.         }
  164.  
  165.         /* Now search for anything making an aspect to the current planet. */
  166.  
  167.         for (j = i+1; j <= cObj; j++) if (!ignore[j] && (fProg || FThing(j)))
  168.           for (k = 1; k <= us.nAsp; k++) if (FAcceptAspect(i, k, j)) {
  169.             d1 = cp1.obj[i]; d2 = cp2.obj[i];
  170.             e1 = cp1.obj[j]; e2 = cp2.obj[j];
  171.             if (MinDistance(d1, d2) < MinDistance(e1, e2)) {
  172.               SwapR(&d1, &e1);
  173.               SwapR(&d2, &e2);
  174.             }
  175.  
  176.             /* We are searching each aspect in turn. Let's subtract the  */
  177.             /* size of the aspect from the angular difference, so we can */
  178.             /* then treat it like a conjunction.                         */
  179.  
  180.             if (MinDistance(e1, Mod(d1-rAspAngle[k])) <
  181.                 MinDistance(e2, Mod(d2+rAspAngle[k]))) {
  182.               e1 = Mod(e1+rAspAngle[k]);
  183.               e2 = Mod(e2+rAspAngle[k]);
  184.             } else {
  185.               e1 = Mod(e1-rAspAngle[k]);
  186.               e2 = Mod(e2-rAspAngle[k]);
  187.             }
  188.  
  189.             /* Check to see if the aspect actually occurs during our    */
  190.             /* segment, making sure we take into account if one or both */
  191.             /* planets are retrograde or if they cross the Aries point. */
  192.  
  193.             f1 = e1-d1;
  194.             if (RAbs(f1) > rDegHalf)
  195.               f1 -= RSgn(f1)*rDegMax;
  196.             f2 = e2-d2;
  197.             if (RAbs(f2) > rDegHalf)
  198.               f2 -= RSgn(f2)*rDegMax;
  199.             if (MinDistance(Midpoint(d1, d2), Midpoint(e1, e2)) < rDegQuad &&
  200.               RSgn(f1) != RSgn(f2)) {
  201.               source[occurcount] = i;
  202.               aspect[occurcount] = k;
  203.               dest[occurcount] = j;
  204.  
  205.               /* Horray! The aspect occurs sometime during the interval.   */
  206.               /* Now we just have to solve an equation in two variables to */
  207.               /* find out where the "lines" cross, i.e. the aspect's time. */
  208.  
  209.               f1 = d2-d1;
  210.               if (RAbs(f1) > rDegHalf)
  211.                 f1 -= RSgn(f1)*rDegMax;
  212.               f2 = e2-e1;
  213.               if (RAbs(f2) > rDegHalf)
  214.                 f2 -= RSgn(f2)*rDegMax;
  215.               g = (RAbs(d1-e1) > rDegHalf ?
  216.                 (d1-e1)-RSgn(d1-e1)*rDegMax : d1-e1)/(f2-f1);
  217.               time[occurcount] = g*divsiz + (real)(div-1)*divsiz;
  218.               sign1[occurcount] = (int)(Mod(cp1.obj[i]+
  219.                 RSgn(cp2.obj[i]-cp1.obj[i])*
  220.                 (RAbs(cp2.obj[i]-cp1.obj[i]) > rDegHalf ? -1 : 1)*
  221.                 RAbs(g)*MinDistance(cp1.obj[i], cp2.obj[i]))/30.0)+1;
  222.               sign2[occurcount] = (int)(Mod(cp1.obj[j]+
  223.                 RSgn(cp2.obj[j]-cp1.obj[j])*
  224.                 (RAbs(cp2.obj[j]-cp1.obj[j]) > rDegHalf ? -1 : 1)*
  225.                 RAbs(g)*MinDistance(cp1.obj[j], cp2.obj[j]))/30.0)+1;
  226.               occurcount++;
  227.             }
  228.           }
  229.       }
  230.     }
  231.  
  232.     /* After all the aspects, etc, in the day have been located, sort   */
  233.     /* them by time at which they occur, so we can print them in order. */
  234.  
  235.     for (i = 1; i < occurcount; i++) {
  236.       j = i-1;
  237.       while (j >= 0 && time[j] > time[j+1]) {
  238.         SwapN(source[j], source[j+1]);
  239.         SwapN(aspect[j], aspect[j+1]);
  240.         SwapN(dest[j], dest[j+1]);
  241.         SwapR(&time[j], &time[j+1]);
  242.         SwapN(sign1[j], sign1[j+1]); SwapN(sign2[j], sign2[j+1]);
  243.         j--;
  244.       }
  245.     }
  246.  
  247.     /* Finally, loop through and display each aspect and when it occurs. */
  248.  
  249.     for (i = 0; i < occurcount; i++) {
  250.       s1 = (int)time[i]/60;
  251.       s2 = (int)time[i]-s1*60;
  252.       j = Day2;
  253.       if (fYear || fProg) {
  254.         l = Mon2;
  255.         while (j > (k = DayInMonth(l, yea0))) {
  256.           j -= k;
  257.           l++;
  258.         }
  259.       }
  260.       SetCI(ciSave, fYear || fProg ? l : Mon, j, yea0,
  261.         DegToDec(time[i] / 60.0), Dst, Zon, Lon, Lat);
  262.       k = DayOfWeek(fYear || fProg ? l : Mon, j, yea0);
  263.       AnsiColor(kRainbowA[k + 1]);
  264.       sprintf(sz, "(%c%c%c) ", chDay3(k)); PrintSz(sz);
  265.       AnsiColor(kDefault);
  266.       sprintf(sz, "%s %s ",
  267.         SzDate(fYear || fProg ? l : Mon, j, yea0, fFalse),
  268.         SzTime(s1, s2)); PrintSz(sz);
  269.       PrintAspect(source[i], sign1[i],
  270.         (int)RSgn(cp1.dir[source[i]])+(int)RSgn(cp2.dir[source[i]]),
  271.         aspect[i], dest[i], sign2[i],
  272.         (int)RSgn(cp1.dir[dest[i]])+(int)RSgn(cp2.dir[dest[i]]),
  273.         (char)(fProg ? 'e' : 'd'));
  274.       PrintInDay(source[i], aspect[i], dest[i]);
  275.     }
  276.   }
  277.   }
  278.  
  279.   /* Recompute original chart placements as we've overwritten them. */
  280.  
  281.   ciCore = ciMain; ciTwin = ciT;
  282.   CastChart(fTrue);
  283. }
  284.  
  285.  
  286. /* Search through a month, year, or years, and print out the times of exact */
  287. /* transits where planets in the time frame make aspect to the planets in   */
  288. /* some other chart, as specified with the -t switch. To do this, we cast   */
  289. /* charts for the start and end of each month, or within a month, and do an */
  290. /* equation check for aspects to the other base chart during the interval.  */
  291.  
  292. void ChartTransitSearch(fProg)
  293. bool fProg;
  294. {
  295.   real planet3[objMax], house3[cSign+1], ret3[objMax], time[MAXINDAY];
  296.   char sz[cchSzDef];
  297.   int source[MAXINDAY], aspect[MAXINDAY], dest[MAXINDAY], sign[MAXINDAY],
  298.     isret[MAXINDAY], M1, M2, Y1, Y2, occurcount, div, i, j, k, s1, s2, s3;
  299.   real divsiz, daysiz, d, e1, e2, f1, f2;
  300.   CI ciT;
  301.  
  302.   ciT = ciTwin;
  303.   for (i = 1; i <= cSign; i++)
  304.     house3[i] = house[i];
  305.   for (i = 1; i <= cObj; i++) {
  306.     planet3[i] = planet[i];
  307.     ret3[i] = ret[i];
  308.   }
  309.  
  310.   /* Hacks: Searching month number zero means to search the whole year    */
  311.   /* instead, month by month. Searching a negative month means to search  */
  312.   /* multiple years, with the span of the year stored in the "day" field. */
  313.  
  314.   Y1 = Y2 = Yea2;
  315.   M1 = M2 = Mon2;
  316.   if (Mon2 < 1) {
  317.     M1 = 1; M2 = 12;
  318.     if (Mon2 < 0) {
  319.       if (Day2 < 1) {
  320.         Y1 = Yea2 + Day2 + 1; Y2 = Yea2;
  321.       } else {
  322.         Y1 = Yea2; Y2 = Yea2 + Day2 - 1;
  323.       }
  324.     }
  325.   }
  326.  
  327.   /* Start searching the year or years in question for any transits. */
  328.  
  329.   for (Yea2 = Y1; Yea2 <= Y2; Yea2++)
  330.  
  331.   /* Start searching the month or months in question for any transits. */
  332.  
  333.   for (Mon2 = M1; Mon2 <= M2; Mon2++) {
  334.     daysiz = (real)DayInMonth(Mon2, Yea2)*24.0*60.0;
  335.     divsiz = daysiz / (real)us.nDivision;
  336.  
  337.     /* Cast chart for beginning of month and store it for future use. */
  338.  
  339.     SetCI(ciCore, Mon2, 1, Yea2, 0.0, Dst2, Zon2, Lon2, Lat2);
  340.     if (us.fProgress = fProg) {
  341.       is.JDp = MdytszToJulian(MM, DD, YY, 0.0, Dst2, Zon2);
  342.       ciCore = ciMain;
  343.     }
  344.     for (i = 1; i <= oNorm; i++)
  345.       SwapN(ignore[i], ignore2[i]);
  346.     CastChart(fTrue);
  347.     for (i = 1; i <= oNorm; i++)
  348.       SwapN(ignore[i], ignore2[i]);
  349.     for (i = 1; i <= cSign; i++)
  350.       cp2.cusp[i] = house[i];
  351.     for (i = 1; i <= oCore; i++) {
  352.       cp2.obj[i] = planet[i];
  353.       cp2.dir[i] = ret[i];
  354.     }
  355.  
  356.     /* Divide our month into segments and then search each segment in turn. */
  357.  
  358.     for (div = 1; div <= us.nDivision; div++) {
  359.       occurcount = 0;
  360.  
  361.       /* Cast the chart for the ending time of the present segment, and */
  362.       /* copy the start time chart from the previous end time chart.    */
  363.  
  364.       d = 1.0 + (daysiz/24.0/60.0)*(real)div/(real)us.nDivision;
  365.       SetCI(ciCore, Mon2, (int)d, Yea2,
  366.         DegToDec(RFract(d)*24.0), Dst2, Zon2, Lon2, Lat2);
  367.       if (fProg) {
  368.         is.JDp = MdytszToJulian(MM, DD, YY, 0.0, Dst2, Zon2);
  369.         ciCore = ciMain;
  370.       }
  371.       for (i = 1; i <= oNorm; i++)
  372.         SwapN(ignore[i], ignore2[i]);
  373.       CastChart(fTrue);
  374.       for (i = 1; i <= oNorm; i++)
  375.         SwapN(ignore[i], ignore2[i]);
  376.       for (i = 1; i <= cSign; i++) {
  377.         cp1.cusp[i] = cp2.cusp[i]; cp2.cusp[i] = house[i];
  378.       }
  379.       for (i = 1; i <= oCore; i++) {
  380.         cp1.obj[i] = cp2.obj[i]; cp1.dir[i] = cp2.dir[i];
  381.         cp2.obj[i] = planet[i];  cp2.dir[i] = ret[i];
  382.       }
  383.  
  384.       /* Now search through the present segment for any transits. Note that */
  385.       /* stars can be transited, but they can't make transits themselves.   */
  386.  
  387.       for (i = 1; i <= cObj; i++) if (!ignore[i])
  388.         for (j = 1; j <= oNorm; j++) if (!ignore2[j] && (fProg || FThing(j)))
  389.  
  390.           /* Between each pair of planets, check if they make any aspects. */
  391.  
  392.           for (k = 1; k <= us.nAsp; k++) if (FAcceptAspect(i, k, j)) {
  393.             d = planet3[i]; e1 = cp1.obj[j]; e2 = cp2.obj[j];
  394.             if (MinDistance(e1, Mod(d-rAspAngle[k])) <
  395.                 MinDistance(e2, Mod(d+rAspAngle[k]))) {
  396.               e1 = Mod(e1+rAspAngle[k]);
  397.               e2 = Mod(e2+rAspAngle[k]);
  398.             } else {
  399.               e1 = Mod(e1-rAspAngle[k]);
  400.               e2 = Mod(e2-rAspAngle[k]);
  401.             }
  402.  
  403.             /* Check to see if the present aspect actually occurs during the */
  404.             /* segment, making sure we check any Aries point crossings.      */
  405.  
  406.             f1 = e1-d;
  407.             if (RAbs(f1) > rDegHalf)
  408.               f1 -= RSgn(f1)*rDegMax;
  409.             f2 = e2-d;
  410.             if (RAbs(f2) > rDegHalf)
  411.               f2 -= RSgn(f2)*rDegMax;
  412.             if (MinDistance(d, Midpoint(e1, e2)) < rDegQuad &&
  413.               RSgn(f1) != RSgn(f2) && occurcount < MAXINDAY) {
  414.  
  415.               /* Ok, we have found a transit. Now determine the time */
  416.               /* and save this transit in our list to be printed.    */
  417.  
  418.               source[occurcount] = j;
  419.               aspect[occurcount] = k;
  420.               dest[occurcount] = i;
  421.               time[occurcount] = RAbs(f1)/(RAbs(f1)+RAbs(f2))*divsiz +
  422.                 (real)(div-1)*divsiz;
  423.               sign[occurcount] = (int)(Mod(
  424.                 MinDistance(cp1.obj[j], Mod(d-rAspAngle[k])) <
  425.                 MinDistance(cp2.obj[j], Mod(d+rAspAngle[k])) ?
  426.                 d-rAspAngle[k] : d+rAspAngle[k])/30.0)+1;
  427.               isret[occurcount] = (int)RSgn(cp1.dir[j]) +
  428.                 (int)RSgn(cp2.dir[j]);
  429.               occurcount++;
  430.             }
  431.           }
  432.  
  433.       /* After all transits located, sort them by time at which they occur. */
  434.  
  435.       for (i = 1; i < occurcount; i++) {
  436.         j = i-1;
  437.         while (j >= 0 && time[j] > time[j+1]) {
  438.           SwapN(source[j], source[j+1]);
  439.           SwapN(aspect[j], aspect[j+1]);
  440.           SwapN(dest[j], dest[j+1]);
  441.           SwapR(&time[j], &time[j+1]);
  442.           SwapN(sign[j], sign[j+1]);
  443.           SwapN(isret[j], isret[j+1]);
  444.           j--;
  445.         }
  446.       }
  447.  
  448.       /* Now loop through list and display all the transits. */
  449.  
  450.       for (i = 0; i < occurcount; i++) {
  451.         s1 = (_int)time[i]/24/60;
  452.         s3 = (_int)time[i]-s1*24*60;
  453.         s2 = s3/60;
  454.         s3 = s3-s2*60;
  455.         SetCI(ciSave, Mon2, s1+1, Yea2, DegToDec((real)
  456.           ((_int)time[i]-s1*24*60) / 60.0), Dst2, Zon2, Lon2, Lat2);
  457.         sprintf(sz, "%s %s ",
  458.           SzDate(Mon2, s1+1, Yea2, fFalse), SzTime(s2, s3)); PrintSz(sz);
  459.         PrintAspect(source[i], sign[i], isret[i], aspect[i],
  460.           dest[i], SFromZ(planet3[dest[i]]), (int)RSgn(ret3[dest[i]]),
  461.           (char)(fProg ? 'u' : 't'));
  462.  
  463.         /* Check for a Solar, Lunar, or any other return. */
  464.  
  465.         if (aspect[i] == aCon && source[i] == dest[i]) {
  466.           AnsiColor(kWhite);
  467.           sprintf(sz, " (%s Return)", source[i] == oSun ? "Solar" :
  468.             (source[i] == oMoo ? "Lunar" : szObjName[source[i]]));
  469.           PrintSz(sz);
  470.         }
  471.         PrintL();
  472. #ifdef INTERPRET
  473.         if (us.fInterpret)
  474.           InterpretTransit(source[i], aspect[i], dest[i]);
  475. #endif
  476.         AnsiColor(kDefault);
  477.       }
  478.     }
  479.   }
  480.  
  481.   /* Recompute original chart placements as we've overwritten them. */
  482.  
  483.   ciCore = ciMain; ciTwin = ciT;
  484.   CastChart(fTrue);
  485. }
  486.  
  487.  
  488. /* Display a list of planetary rising times relative to the local horizon */
  489. /* for the day indicated in the chart information, as specified with the  */
  490. /* -Zd switch. For the day, the time each planet rises (transits horizon  */
  491. /* in East half of sky), sets (transits horizon in West), reaches its     */
  492. /* zenith point (transits meridian in South half of sky), and nadirs      */
  493. /* transits meridian in North), is displayed.                             */
  494.  
  495. void ChartInDayHorizon()
  496. {
  497.   char sz[cchSzDef];
  498.   int source[MAXINDAY], type[MAXINDAY], sign[MAXINDAY],
  499.     fRet[MAXINDAY], occurcount, division, div, s1, s2, i, j, fT;
  500.   real time[MAXINDAY], rgalt1[objMax], rgalt2[objMax], azialt[MAXINDAY],
  501.     azi1, azi2, alt1, alt2, lon, lat, mc1, mc2, xA, yA, xV, yV, d, k;
  502.   CI ciT;
  503.  
  504.   fT = us.fSidereal; us.fSidereal = fFalse;
  505.   lon = RFromD(Mod(Lon)); lat = RFromD(Lat);
  506.   division = us.nDivision * 4;
  507.   occurcount = 0;
  508.  
  509.   ciT = ciTwin; ciCore = ciMain; ciCore.tim = 0.0;
  510.   CastChart(fTrue);
  511.   mc2 = RFromD(planet[oMC]); k = RFromD(planetalt[oMC]);
  512.   EclToEqu(&mc2, &k);
  513.   for (i = 1; i <= cSign; i++) {
  514.     cp2.cusp[i] = house[i];
  515.     cp2.house[i] = inhouse[i];
  516.   }
  517.   for (i = 1; i <= cObj; i++) {
  518.     cp2.obj[i] = planet[i];
  519.     rgalt2[i] = planetalt[i];
  520.     cp2.dir[i] = ret[i];
  521.   }
  522.  
  523.   /* Loop thorough the day, dividing it into a certain number of segments. */
  524.   /* For each segment we get the planet positions at its endpoints.        */
  525.  
  526.   for (div = 1; div <= division; div++) {
  527.     ciCore = ciMain; ciCore.tim = DegToDec(24.0*(real)div/(real)division);
  528.     CastChart(fTrue);
  529.     mc1 = mc2;
  530.     mc2 = RFromD(planet[oMC]); k = RFromD(planetalt[oMC]);
  531.     EclToEqu(&mc2, &k);
  532.     for (i = 1; i <= cSign; i++) {
  533.       cp1.cusp[i] = cp2.cusp[i]; cp1.house[i] = cp2.house[i];
  534.       cp2.cusp[i] = house[i];    cp2.house[i] = inhouse[i];
  535.     }
  536.     for (i = 1; i <= cObj; i++) {
  537.       cp1.obj[i] = cp2.obj[i]; cp2.obj[i] = planet[i];
  538.       rgalt1[i] = rgalt2[i]; rgalt2[i] = planetalt[i];
  539.       cp1.dir[i] = cp2.dir[i]; cp2.dir[i] = ret[i];
  540.     }
  541.  
  542.     /* For our segment, check to see if each planet during it rises, sets, */
  543.     /* reaches its zenith, or reaches its nadir.                           */
  544.  
  545.     for (i = 1; i <= cObj; i++) if (!ignore[i] && FThing(i)) {
  546.       EclToHorizon(&azi1, &alt1, cp1.obj[i], rgalt1[i], lon, lat, mc1);
  547.       EclToHorizon(&azi2, &alt2, cp2.obj[i], rgalt2[i], lon, lat, mc2);
  548.       j = 0;
  549.  
  550.       /* Check for transits to the horizon. */
  551.       if ((alt1 > 0.0) != (alt2 > 0.0)) {
  552.         d = RAbs(alt1)/(RAbs(alt1)+RAbs(alt2));
  553.         k = Mod(azi1 + d*MinDifference(azi1, azi2));
  554.         j = 1 + 2*(MinDistance(k, rDegHalf) < rDegQuad);
  555.  
  556.       /* Check for transits to the meridian. */
  557.       } else if (RSgn(MinDifference(azi1, rDegQuad)) !=
  558.         RSgn(MinDifference(azi2, rDegQuad))) {
  559.         j = 2 + 2*(MinDistance(azi1, rDegQuad) < rDegQuad);
  560.         d = RAbs(azi1 - (j > 2 ? rDegQuad : 270.0))/MinDistance(azi1, azi2);
  561.         k = alt1 + d*(alt2-alt1);
  562.       }
  563.       if (j && occurcount < MAXINDAY) {
  564.         source[occurcount] = i;
  565.         type[occurcount] = j;
  566.         time[occurcount] = 24.0*((real)(div-1)+d)/(real)division*60.0;
  567.         sign[occurcount] = (int)Mod(cp1.obj[i] +
  568.           d*MinDifference(cp1.obj[i], cp2.obj[i]))/30 + 1;
  569.         fRet[occurcount] = (int)RSgn(cp1.dir[i]) + (int)RSgn(cp2.dir[i]);
  570.         azialt[occurcount] = k;
  571.         ciSave = ciMain;
  572.         ciSave.tim = DegToDec(time[occurcount] / 60.0);
  573.         occurcount++;
  574.       }
  575.     }
  576.   }
  577.  
  578.   /* Sort each event in order of time when it happens during the day. */
  579.  
  580.   for (i = 1; i < occurcount; i++) {
  581.     j = i-1;
  582.     while (j >= 0 && time[j] > time[j+1]) {
  583.       SwapN(source[j], source[j+1]);
  584.       SwapN(type[j], type[j+1]);
  585.       SwapR(&time[j], &time[j+1]);
  586.       SwapN(sign[j], sign[j+1]);
  587.       SwapN(fRet[j], fRet[j+1]);
  588.       SwapR(&azialt[j], &azialt[j+1]);
  589.       j--;
  590.     }
  591.   }
  592.  
  593.   /* Finally display the list showing each event and when it occurs. */
  594.  
  595.   for (i = 0; i < occurcount; i++) {
  596.     ciSave = ciMain;
  597.     ciSave.tim = DegToDec(time[i] / 60.0);
  598.     j = DayOfWeek(Mon, Day, Yea);
  599.     AnsiColor(kRainbowA[j + 1]);
  600.     sprintf(sz, "(%c%c%c) ", chDay3(j)); PrintSz(sz);
  601.     AnsiColor(kDefault);
  602.     s1 = (int)time[i]/60;
  603.     s2 = (int)time[i]-s1*60;
  604.     sprintf(sz, "%s %s ", SzDate(Mon, Day, Yea, fFalse), SzTime(s1, s2));
  605.     PrintSz(sz);
  606.     AnsiColor(kObjA[source[i]]);
  607.     sprintf(sz, "%7.7s ", szObjName[source[i]]); PrintSz(sz);
  608.     AnsiColor(kSignA(sign[i]));
  609.     sprintf(sz, "%c%c%c%c%c ",
  610.       fRet[i] > 0 ? '(' : (fRet[i] < 0 ? '[' : '<'), chSig3(sign[i]),
  611.       fRet[i] > 0 ? ')' : (fRet[i] < 0 ? ']' : '>')); PrintSz(sz);
  612.     AnsiColor(kElemA[type[i]-1]);
  613.     if (type[i] == 1)
  614.       PrintSz("rises  ");
  615.     else if (type[i] == 2)
  616.       PrintSz("zeniths");
  617.     else if (type[i] == 3)
  618.       PrintSz("sets   ");
  619.     else
  620.       PrintSz("nadirs ");
  621.     AnsiColor(kDefault);
  622.     PrintSz(" at ");
  623.     if (type[i] & 1) {
  624.       j = (int)(RFract(azialt[i])*60.0);
  625.       sprintf(sz, "%3d%c%02d'", (int)azialt[i], chDeg1, j); PrintSz(sz);
  626.  
  627.       /* For rising and setting events, we'll also display a direction  */
  628.       /* vector to make the 360 degree azimuth value thought of easier. */
  629.  
  630.       xA = RCosD(azialt[i]); yA = RSinD(azialt[i]);
  631.       if (RAbs(xA) < RAbs(yA)) {
  632.         xV = RAbs(xA / yA); yV = 1.0;
  633.       } else {
  634.         yV = RAbs(yA / xA); xV = 1.0;
  635.       }
  636.       sprintf(sz, " (%.2f%c %.2f%c)",
  637.         yV, yA < 0.0 ? 's' : 'n', xV, xA > 0.0 ? 'e' : 'w'); PrintSz(sz);
  638.     } else
  639.       PrintAltitude(azialt[i]);
  640.     PrintL();
  641.   }
  642.  
  643.   /* Recompute original chart placements as we've overwritten them. */
  644.  
  645.   ciCore = ciMain; ciTwin = ciT;
  646.   us.fSidereal = fT;
  647.   CastChart(fTrue);
  648. }
  649.  
  650.  
  651. /* Print out an ephemeris - the positions of the planets (at the time in the */
  652. /* current chart) each day during a specified month, as done with the -E     */
  653. /* switch. Display the ephemeris for the whole year if -Ey is in effect.     */
  654.  
  655. void ChartEphemeris()
  656. {
  657.   char sz[cchSzDef];
  658.   int yea, yea1, yea2, mon, mon1, mon2, daysiz, i, j, s, d, m;
  659.  
  660.   /* If -Ey is in effect, then loop through all months in the whole year. */
  661.  
  662.   if (us.nEphemYears) {
  663.     yea1 = Yea; yea2 = yea1 + us.nEphemYears - 1; mon1 = 1; mon2 = 12;
  664.   } else {
  665.     yea1 = yea2 = Yea; mon1 = mon2 = Mon;
  666.   }
  667.  
  668.   /* Loop through the year or years in question. */
  669.  
  670.   for (yea = yea1; yea <= yea2; yea++)
  671.  
  672.   /* Loop through the month or months in question, printing each ephemeris. */
  673.  
  674.   for (mon = mon1; mon <= mon2; mon++) {
  675.     daysiz = DayInMonth(mon, yea);
  676.     PrintSz(us.fEuroDate ? "Dy/Mo/Yr" : "Mo/Dy/Yr");
  677.     for (j = 1; j <= cObj; j++) {
  678.       if (!ignore[j] && FThing(j)) {
  679.         sprintf(sz, "  %s%c%c%c%c", is.fSeconds ? "  " : "", chObj3(j),
  680.           szObjName[j][3] != 0 ? szObjName[j][3] : ' '); PrintSz(sz);
  681.         PrintTab(' ', us.fParallel ? 2 + is.fSeconds : 1 + 3*is.fSeconds);
  682.       }
  683.     }
  684.     PrintL();
  685.     for (i = 1; i <= daysiz; i = AddDay(mon, i, yea, 1)) {
  686.  
  687.       /* Loop through each day in the month, casting a chart for that day. */
  688.  
  689.       SetCI(ciCore, mon, i, yea, Tim, Dst, Zon, Lon, Lat);
  690.       CastChart(fTrue);
  691.       PrintSz(SzDate(mon, i, yea, -1));
  692.       PrintCh(' ');
  693.       for (j = 1; j <= cObj; j++)
  694.         if (!ignore[j] && FThing(j)) {
  695.           if (!us.fParallel) {
  696.             if (is.fSeconds)
  697.               PrintZodiac(planet[j]);
  698.             else {
  699.               AnsiColor(kObjA[j]);
  700.               s = SFromZ(planet[j]);
  701.               d = (int)planet[j] - (s-1)*30;
  702.               m = (int)(RFract(planet[j])*60.0);
  703.               sprintf(sz, "%2d%s%02d", d, szSignAbbrev[s], m); PrintSz(sz);
  704.             }
  705.           } else {
  706.             AnsiColor(kObjA[j]);
  707.             PrintAltitude(planetalt[j]);
  708.           }
  709.           PrintCh((char)(ret[j] >= 0.0 ? ' ' : '.'));
  710.         }
  711.       PrintL();
  712.       AnsiColor(kDefault);
  713.     }
  714.     if (mon < mon2 || yea < yea2)
  715.       PrintL();
  716.   }
  717. }
  718.  
  719. /* charts3.c */
  720.